Particle-hole symmetry and interaction effects in the Kane-Mele-Hubbard model 
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We prove that the Kane-Mele-Hubbard model with purely imaginary next-nearest-neighbor hop- 
pings has a particle-hole symmetry at half-filling. Such a symmetry has interesting consequences 
including the absence of charge and spin currents along open edges, and the absence of the sign 
problem in the determinant quantum Monte- Carlo simulations. Consequentially, the interplay be- 
tween band topology and strong correlations can be studied at high numeric precisions. The process 
that the topological band insulator evolves into the antiferromagnetic Mott insulator as increasing 
interaction strength is studied by calculating both the bulk and edge electronic properties. In agree- 
ment with previous theory analyses, the numeric simulations show that the Kane-Mele-Hubbard 
model exhibits three phases as increasing correlation effects: the topological band insulating phase 
with stable helical edges, the bulk paramagnetic phase with unstable edges, and the bulk antiferro- 
magnetic phase. 
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I. INTRODUCTION 

The precise quantization of the Hall conductance in 
the integer quantum Hall states is protected by the 
non-trivial topology of band structures. This topologi- 
cal property is characterized by the Thouless-Kohmoto- 
Nightingale-den Nijs (TKNN) number, or the Chern 
number p], Q, which takes non-zero values only when 
time-reversal symmetry is broken. In recent years, 
tremendous progress has been achieved in a new class of 
topologically non-trivial band insulators in the presence 
of time-reversal symmetry, which are termed as topologi- 
cal insulators (3l— Tl3l] . Topological insulators exist in both 
two (2D) and three dimensions (3D), which are charac- 
terized by the Z2 topological index. These topological 
states have robust gapless helical edge modes with odd 
number of channels in 2D 0JII,[l5], and odd number of 
surface Dirac cones in 3D [lll-[l3j . Topological insulators 
have been experimentally observed in 2D quantum wells 
through transport measurements [l6[ ? and also in 3D sys- 
tems of Bi^Sbi-^, Bi2Te3, Bi2Se3, and Sb2Te3 through 
the angle-resolved photo-emission spectroscopy [17H2(J, 
and the absence of backscattering in the scanning tun- 
neling spectroscopy [2~0-[23|]. 

Interaction effects in topological insulators remain an 
open question. Due to their gapped nature, topologi- 
cal insulators remain stable against weak interactions. 
However, strong interactions may change their topolog- 
ical properties. For 2D topological insulators, it has 
been found that the two-particle correlated backscatter- 
ing, which is an interaction effect and is allowed in the 
time-reversal invariant Hamiltonian, can gap out the he- 
lical edge states by spontaneous developing magnetic or- 
dering under strong repulsive interactions |l4l [IE). In 
this case, time-reversal symmetry is spontaneously bro- 
ken along edges, although the bulk remains paramag- 
netic. At mean-field level, interaction effects can desta- 
bilize the quantum anomalous Hall state of the Haldane- 



Hubbard model [24| and the 2D topological insulating 
state of the Kane-Mele-Hubbard (KMH) model [H[ by 
developing long-range charge density wave and antifer- 
romagnetic orders, respectively [25j. Interactions can 
also change the topologically trivial band structures into 
non-trivial ones at mean-field level by developing bulk 
order parameters [261429^ Due to the difficulty of ana- 
lytic studies on strong correlation physics, exact results 
from numeric simulations are desirable. Recently, an 
exact diagonalization has been carried on the spinless 
Haldane-Hubbard model [30]. A first order phase tran- 
sition between quantum anomalous Hall insulating state 
and topologically trivial Mott-insulating state is found. 

Quantum Monte- Carlo (QMC) simulations play an im- 
portant role in studying strongly correlated systems [3ll — 
134 ]. A major obstacle to apply the QMC to fermion 
systems is the notorious sign problem. In the partic- 
ular method of the determinant QMC, the 4- fermion 
interaction terms are decoupled through the Hubbard- 
Stratonovich (HS) transformation and fermions are able 
to be integrated out. The resultant fermion determinant, 
generally speaking, is not positive-definite, which is the 
origin of the notorious sign problem. This problem pre- 
vents QMC simulations to achieve a good numerical pre- 
cision at low temperatures and large sample sizes. Never- 
theless, in a number of interacting models, the sign prob- 
lem disappears. As presented in Ref. [35|, these models 
include the negative- U Hubbard model, the positive- U 
Hubbard model at half-filling and in bipartite lattices, 
and a class of models whose interactions can be decom- 
posed in a time-reversal invariant way. 

We find that the Kane-Mele model augmented by the 
Hubbard interaction with purely imaginary next-nearest- 
neighbor hoppings has a particle-hole symmetry. Such a 
symmetry has interesting consequences such as the ab- 
sence of edge charge and spin currents, which shows the 
edge currents are not a reliable criterion for topological 
properties. More importantly, the particle- hole symme- 
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try ensures the absence of the sign problem in the quan- 
tum Monte-Carlo simulations. This provides a wonderful 
opportunity to study interaction effects in topological in- 
sulating systems. 

In this article, we perform a determinant QMC study 
on the stability of the topological insulating state of the 
KMH model with the strong Hubbard interaction U . An- 
tiferromagnetic long-range-order has been found at large 
values of U. Consequently, the quantum phase diagram 
of the KMH model can be classified into paramagnetic 
bulk insulating phases and antiferromagnetic Mott in- 
sulating phases. When further consider the stability of 
helical edges with infinitesimal two-particle backscatter- 
ing, which is not contained in KMH model but gen- 
erally allowed by time-reversal symmetry, the param- 
agnetic bulk insulating phase can be divided into two 
regimes according to their edge state Luttinger param- 
eters [14]. The topological band insulator with stable 
helical edges are stable in the weak interaction regime, 
while the helical edges become unstable by two-particle 
correlated backscattering at the intermediate interaction 
regime. We have also studied the nature of spin-liquid 
phase in the pure Hubbard model with A = 0, showing 
that it is neither a spontaneous Haldane type quantum 
anomalous Hall insulator, nor, a Kane-Mele type quan- 
tum spin Hall insulator. 

This article is organized as follows. In Section [IIJ we 
prove the absence of the sign problem in the KMH model 
under certain conditions. In Section HIH we present the 
simulations on the developing of antiferromagnetic long- 
range orders in the bulk. In Section [IVJ the edge prop- 
erties are studied including both the edge single particle 
excitations and the edge spin correlations. In Section W\ 
we present the simulation of the charge and spin current 
orders in the pure Hubbard model in the honeycomb lat- 
tice. Conclusions are given in Section fVTl 



II. GENERAL PROPERTIES OF THE KMH 
MODEL 

The Kane-Mele model is a straightforward generaliza- 
tion of the Haldane model in the honeycomb lattice Q 
defined as 

(i,3),<r ((i t i'))a,P 
~ 4'a a z,apCip} ~ /i^cJ^C^, (1) 

where t is the nearest-neighbor (NN) hopping integral as 
scaled to 1 below; A is the next-nearest-neighbor (NNN) 
spin-orbit hopping integral; \i is the chemical potential. 
In the general case of the Kane-Mele model, the NNN 
hopping for the spin-^ (I) electrons are complex-valued 
and complex-conjugate to each other. As a special case, 
the NNN hopping in Eq. [1] is purely imaginary. The 



Hubbard interaction is defined as usual 

Hint = uY,[^-\][nn-\]. (2) 

i 

In this section, we will present the symmetry properties 
of Eq. [1] and Eq. [2j and prove the absence of the sign 
problem in the determinant QMC. 

A. Particle-hole symmetry 

Eq. [D and Eq. [2] has the particle-hole symmetry at 
H = as explained below. We define the transformation 
as usual 




dia = (-l)'cL — ► 4r = ("l)^- (3) 



Under this transformation, a Hermitian fermion bilinear 
operator connecting two sites belonging to two different 
sublattices transforms as 

c \(jKijCj a + Cj a (Kij) Ci<j > 

dUKijYd^+dlKijd^, (4) 

while that connecting two different sites in the same sub- 
lattice transforms as 

cl a KurCi' a +cl, a (Kiit)*Cit„ — > 
- dl(K iif yd ifa -4 a K iif d ifa . (5) 

The onsite particle density transforms as 

where no summation over spin-index is assumed in Eq. 
[6l Clearly in Eq. [TJ the NN-hopping is real and the 
NNN-hopping is purely imaginary, thus its band struc- 
ture is invariant at \i = 0. Eq. [2] is obviously invariant. 
The particle-hole symmetry also implies that \i = cor- 
responds to half-filling. 

B. Absence of the charge and spin currents 

An important conclusion based on the particle-hole 
symmetry is that both charge and spin currents vanish 
on all the bonds for the KMH model of Eq. [Hand Eq. [2] 
at \i = 0. This result applies to arbitrary boundary con- 
ditions with broken bonds but with the homogeneous on- 
site potential which maintains the particle-hole symme- 
try on each site. The proof is straightforward. Through 
the continuity equation, the current operators of each 
spin component along the NN and NNN bonds are de- 
fined as 

jNN t c . _J c . \ 
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respectively, where no summation over spin-index is as- 
sumed. Both J NN and J NNN are odd under the particle- 
hole transformation, thus they vanish even with the open- 
boundary condition. By the same reasoning, the charge 
current also vanishes in the Haldane-Hubbard model with 
the purely imaginary NNN-hoppings and the particle- 
hole symmetric charge interactions of 



HNN,int = ^2 ViA n i 
ij 



(8) 



This result shows that edge charge and spin currents 
are not good criteria for quantum anomalous Hall and 
topological insulators. In order to have a better under- 
standing on this counter-intuitive result, we have con- 
sidered the simplest non-interacting Haldane model with 
the purely imaginary NNN hoppings by diagonalization. 
There are indeed gapless one-dimensional single particle 
chiral edge modes clearly seen from the spectra as com- 
monly presented in literatures. Clearly this branch of 
edge mode contributes to edge currents. However, we 
find that the continuous bulk spectra also contribute to 
edge currents. Perfect cancellation occurs which results 
in zero current on each bond, including each edge bond, 
although we know for sure that the band structure is 
topologically non-trivial. For interacting models, there 
are no well-defined single particle states. We cannot sep- 
arate the edge and bulk contributions anymore. Nev- 
ertheless, we expect that current correlation functions 
should exhibit difference between topological insulators 
and trivial insulators. 

Another conclusion inferred from the particle-hole 
symmetry is that the average particle density for each 
spin component on each site is strictly \ even when the 
translational symmetry is broken. For example, it ap- 
plies to any disordered pattern of the hopping integrals, 
as long as the NN hoppings are real and the NNN hop- 
pings are purely imaginary. 

Edge currents do appear if the particle-hole symme- 
try is broken. For example, for the non-interacting Hal- 
dane model with generally complex-valued NNN hop- 
pings, edge currents appear along open boundaries. So 
far we only consider the sharp edges of broken bonds but 
with homogeneous on-site potential. For edges with the 
confining single particle potential, the particle-hole sym- 
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metry is broken which also results in edge currents. In 
particular, for a weak linear external potential, the linear 
response should still give rise to quantized Hall conduc- 
tance in the insulating region. 



C. Absence of the QMC sign problem 

The Hubbard model on the honeycomb lattice, which 
corresponding to the case of A = of Eq. [I] and Eq. [2] has 
been recently simulated at half-filling [36|. As Hubbard 
U increases from zero to a moderate value and then the 
strong coupling regime, the ground state emerges from a 
semi- metal phase, to a new spin-liquid phase and then the 
antiferromagnetic insulating phase. Below we will prove 
that the sign problem still vanishes with nonzero values 
of A. The absence of the sign problem can be proved for 
both the finite temperature and the zero temperature al- 
gorithms for the determinant QMC. In this subsection, 
we prove this property for the finite temperature method 
for simplicity, and leave the more lengthy proof for the 
zero temperature algorithm in Appendix [XJ We empha- 
size that the simulations presented in this article are done 
at the zero temperature. 

Just as the Ref.([36[) does, we employ a discrete HS 
transformation which respects the SU(2) symmetry for 
every fixed HS field configuration by decoupling in the 
density channel. We rewrite the Hubbard interaction and 
decompose it in the density channel by using imaginary 
numbers as 



-AU(n t +n l -l) 2 /2 



^2 7 i (/)e i ^ (z) V^ r ¥( n t+^- 1 ) 

l=±l,±2 

0(Ar 4 ). (9) 



where the discretized HS fields take values of 7(±1) = 
1 + x/6/3, 7 (±2) = 1 - V6/3; rj(±l) = ±^2(3-^6), 

and rj(±2) = ±^2(3 + ^). 

For the convenience of presentation, we prove the ab- 
sence of the sign problem in the finite temperature for- 
malism with p = 1/T. The proof for the zero tempera- 
ture projector algorithm is similar. The partition func- 
tion at half-filling reads 



{1} \ p=M 




^E ilj cl ; 4c j ; e iVAr[//2E i % lP (0(c! ; c i ;-i) 



(10) 



where sums over all the configurations of the dis- crete HS fields 7]i iP (l) and 7i, p (Z); i and p are indices of dis- 
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cretized grids along the spatial and temporal directions, 
respectively; Tr takes the trace of the fermion space; At 
is the discretized time slice which is set to 0.05 in the sim- 
ulations in this article; MAr equals the imaginary time 
p. By using the particle-hole transformation defined in 
Eq. [3j we show that the onsite particle density trans- 
forms according to Eq. [6j the NN-hopping matrix kernel 
transforms according to Eq. HJ the NNN-hopping matrix 
kernel transforms according to Eq. [51 

When the following two conditions are satisfied, 
the fermion determinants of two spin components are 
complex-conjugate to each other, thus the product of 
them is positive-definite: 

K% = (Kfi)* = Kfj for NN-hopping; 

K% = -(KJ i )* = -Kf j for NNN-hopping. (11) 

Apparently, Eq. [1] and Eq. [2] satisfy these conditions, 
and thus are sign problem free. 

Please note that the KMH mode is sign problem free 
only when the NNN-hopping is purely imaginary. Gen- 
erally speaking, the interacting model without the sign 
problem can have complex-valued hoppings with oppo- 
site signs, which still gives rise to opposite Chern num- 
bers for the band structures of spin-^ and ^, respectively. 
However, they are not related by time-reversal symmetry 
anymore. 
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FIG. 1: (Color online) The comparison between the antifer- 
romagnetic structure factors S^f along the z-axis and S^f 
in the xy-p\ane at A = 0.1 for the size oiN — 2xLxL with 
L — 6. The easy-plane feature is clear. 



model at zero temperature by using the projective method 
[371 ]. We perform measurements from 10 different random 
number series and each independent measurement has 
500 sample sweeps after warming up, the discrete imagi- 
nary time step Ar is set to be 0.05. In this section, we use 
periodic boundary conditions for bulk properties calcu- 
lation, e.g., the bulk antiferromagnetic structure factor. 



III. THE QMC STUDY ON THE BULK 
PROPERTIES OF THE KMH-MODEL 

The Hubbard model in the honeycomb lattice, which 
corresponds the case of A = in Eq. [1] and Eq. [2j has 
been simulated in Ref. [36]. When U increases from 
zero, the single particle charge gap appears at U = 3.7, 
while the antiferromagnetic long-rang order emerges at 
U = 4.3. The mismatch reveals an exotic spin liquid 
phase in between. When the intrinsic spin-orbit coupling, 
i.e., the NNN hopping term in Eq. [TJ enters, the model 
describes the topological band insulator. It already has 
a band gap even at U = 0. As increasing U, the anti- 
ferromagnetic structure factor is still a good quantity to 
tell when the magnetic long range order appears. How- 
ever, the bulk gap is no longer an appropriate quantity 
to judge a possible transition from the topological band 
insulator to a antiferromagnetic Mott- insulator. Here we 
use the local single particle gap on edge sites as an indica- 
tor of the stability of edge states and topological proper- 
ties. We also study the edge effects to antiferromagnetic 
correlations. In this section, we will simulate the bulk 
antiferromagnetic structure factor, and leave the study 
of edge properties in Section IIVI 



A. Sampling parameters of our simulations 

Based on the above proof of the absence of the sign 
problem, we perform the QMC simulation for the KMH 



B. The developing of the bulk antiferromagnetic 
long range order 

The spin-orbit NNN hopping in Eq. 1 breaks the 
SU(2) symmetry but preserves the conservation of S z . 
As a result, the antiferromagnetic correlation of S z 
should be different from those of S x and S y . In the large 
[/-limit, the NNN hopping generates an anisotropic ex- 
change as 



T/ / QX QX j_ qV qV qZqZ\ 



(12) 



with J' = 4A 2 //7, which is ferromagnetic in the xy-plane 
and antiferromagnetic along the z-direction (25[. As the 
combined effect from the NNN anisotropic exchange and 
NN isotropic antiferromagnetic exchange, the magnetic 
exchange along the z-axis is frustrated while those along 
x and y- axes are not. Thus the Neel ordering favors the 
easy x?/-plane. 

Our QMC simulations have confirmed this picture. 
The antiferromagnetic structure factor along the in- 
direction (xx-AFSF) and the ^-direction (zz-AFSF) are 
defined as 



STf = jj(G\ 



E(- 1 ) i ^ z 



\G), 



\G), 



(13) 




FIG. 2: (Color online) The finite-size scaling of the xx- 
antiferromagnetic structure factors calculated at A = 0.1 for 
the sizes of N = 2 x L x L (L = 3, 6, 9 and 12), and the 
different values of U indicated in the inset. Finite values of 
Saf/N in the thermodynamic limit appear at U > U c with 
Uc ~ 4.9. 




FIG. 3: (Color online) The QMC simulation of the phase dia- 
gram of the KMH model. The antiferromagnetically long- 
range ordered phase appears at strong correlation regime. 
The paramagnetic phase is divided into two regimes: topolog- 
ical band insulator (TBI) with stable helical edges, and bulk 
paramagnetic phase with unstable edges (see further discus- 
sions in Sect. IIV C|) . The two critical values of U at A = 
are from Ref. [36|] by Meng et al, which are also confirmed 
in our QMC simulations. 



where (G\..\G) means average over the ground state; N = 
2 x L x L is the number of sites; L is the size; (— )* takes 
the values of ±1 for the A and 5-sublattices, respectively. 
The comparison between Sjf F and S y ^ F is plotted in Fig. 
[TJ which clearly shows the easy-plane feature. 

Below we will use the xx-AFSF to describe the anti- 
ferromagnetic properties, and perform the simulation at 
A = 0.1 with different values of U and sample sizes of 
L = 3, 6, 9, 12. The extrapolation to the thermodynamic 
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zigzag edge 
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FIG. 4: (Color online) The KMH model in the lattice with 
the zig-zag edges. The boundary conditions are periodical and 
open along the x and ^-directions, respectively. The NNN- 
bonds between two closest tips on the zig-zag edges are re- 
moved. 



limit for different Hubbard U is plotted in of Fig. El It 
can be seen that the magnetic long range order emerges 
at U c = 4.9 ± 0.1 for A = 0.1. In Fig. [3] we present the 
QMC simulation on the magnetic phase diagram of the 
KMH model for in the parameter space of (U, A). The 
phase boundary separating the AF long-range-ordered 
phase and non-magnetic phases are marked for various 
values of A. The spin-orbit coupling opens the band gap 
at the order of A, thus the interaction effect U becomes 
important only when U is larger than A. As a result, the 
critical value of U c for the onset of the AF phase increases 
with A. 

The phase diagram Fig. [3] exhibits a large regime 
of non-magnetic insulating state outside the AF phase 
at A 7^ 0. At small values of ?7, it should be the Z2 
topological band insulating phase which is stable against 
weak interactions. As increasing U, it enters the AF 
Mott insulating phase at a critical line of U c . In an up- 
dated version of Ref. [40] , it is found that the spin-liquid 
phase also extends to a small but finite value of A. How- 
ever, the nature of this spin-liquid state remains unclear. 
The bulk paramagnetic regime actually has rich internal 
structures. According to the stability of the helical edge 
states with respect to the two-particle spin-flip backscat- 
tering, this paramagnetic insulating phase is divided into 
two different regimes with the effective edge Luttinger 
parameter K < (>)|, respectively. The analysis is pre- 
sented below in Sect. IIV CI 



IV. THE QMC STUDY OF THE EDGE 
PROPERTIES OF THE KMH MODEL 

We believe that the edge properties is crucial to expose 
the topological aspect of the KMH model. In this sec- 
tion, we will show that the antiferromagnetic correlations 
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FIG. 5: (Color online) The extrapolation of local single par- 
ticle gap for the tip sites on the zig-zag edges of a ribbon 
geometry with the size of 2 x L x n y with n y — 8. In the in- 
set, the logarithms of onsite time-displaced Green's functions 
lnG(i, i;r) of the tip sites is depicted for U = 2. The slopes 
of the long time tails measure the edge excitation gap A e d ge - 
Here A is set to be 0.1 in this calculation. Due to the one- 
dimensional nature of the edge, its densities of states in the 
thermodynamic limit are depleted according to the power-law 
determined by the Luttinger parameter, and thus are difficult 
to be distinguished from finite gaps. 



along the edge become strongly relevant as increasing the 
Hubbard U while the bulk remains paramagnetic. We 
consider the lattice configuration plotted in Fig. 2] with 
the periodical and open boundary conditions along the x 
and ^/-directions, respectively. 



A. The single-particle excitations 

As proved in Sect. [Ill the edge currents, both for 
charge and spin, are always zero due to the particle-hole 
symmetry. We use another quantity, the local single par- 
ticle excitation gap on edge sites, to check whether the 
edges are gapped or gapless. It can be extracted from 
the tail of on-site time displaced Green's function on the 
edge lnG(i,i;r) ~ A e ^ e r, which is defined by 

G(i,i;r) = ±<G| £ cj t (r)c it (0) + c^(r) Cii (0)|G), 



itztip 



(14) 



where \G) is the many-body ground state. The depen- 
dency of In G(i : i;r) with r for the site i on the tip of 
the zig-zag edges are plotted in the inset of Fig. [5j where 
the long tail of In G(z,z;r) shows a linear behavior with 
r and the slope measures the excitation gap. Here the 
lattice has a ribbon geometry with n y zig-zag rows. We 
fix the width of the ribbon n y = 8 and increase its length. 
The extrapolations of the edge excitation gaps with L are 
depicted in Fig. [5] with A fixed at 0.1 and different values 




3 4 5 6 7 
Index of zigzag rows 

FIG. 6: (Color online) The row xx-AFSF defined in Eq. fTBlfor 
each zig-zag rows parallel to the boundary. The parameters 
values are A = 0.1; the sample size 2 x L x L with L = 8; and 
different values of U indicated in the inset. The row- indices 
1 and 8 correspond to the boundary rows, and those of 4 and 
5 corresponds to the central rows. 



ofU < U c . Clearly increasing U significantly reduces the 
weight of the low energy spectra. 

The bosonization analysis of the stability of the helical 
edge states has been performed in Ref. [14], LL5[ . For the 
parameter regime of Fig. [5j the bulk remains paramag- 
netic, or, time-reversal invariant. For the current KMH- 
model, S z is conserved which prohibits the existence of 
the two-particle spin-flip scattering term to open the gap. 
The Luttinger liquid theory of such a helical edge branch, 
i.e., the right and left movers are with opposite spin po- 
larizations, is characterized by only one Luttinger param- 
eter K, which describes the forward scattering between 
these two branches. Due to the helical nature of the edge 
states, the long wavelength charge fluctuations and the 
z-component of the spin fluctuations are not indepen- 
dent but are conjugate to each other. Both of them are 
gapless in the thermodynamic limit, and so does the sin- 
gle particle edge excitations. The onsite imaginary time 
single-particle Green's function decays as l/r a with the 
exponent 



a = K+l/K. 



(15) 



At K <C 1, the low energy density of states does not open 
a full gap but are depleted according to a power-law, and 
thus exhibit a pseudo-gap behavior. The non-zero gap 
values in Fig. [5] may be an artifact of finite size scaling 
and a result of tunneling between two opposite edges. 
A more detailed numerical analysis is needed to further 
clarify the nature of the single particle excitations. 



B. Edge spin structure factors 

We further investigate the edge effects to the antiferro- 
magnetic correlations. We define the antiferromagnetic 
structure form factor for each zig-zag row parallel to the 
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FIG. 7: (Color online) The finite size scaling of the xx-AFSF 
defined in Eq. [16] for the edge row with A = 0.1. The size of 
this ribbon is 2 x L x 4. We emphasize that due to the ID 
nature of the edge and the U(l) spin symmetry, this scaling 
actually shows the power-law correlation rather than the true- 
long-range order. The finite intercepts are mainly due to small 
size effects. 



zig-zag boundary as 



1 Zigzag, AF 



(m) 



±<G|[£(-1)^/|G>, 



(16) 



where m is the index of the zig-zag row; i is the site index 
along the m-th zig-zag line; 2L is the number of sites in 
each row. The xx-AFSF for all the rows are depicted in 
Fig. El 

It is interesting to observe that the AF correlations are 
strongest on edges, and become weaker inside the bulk. 
This effect is most prominent at small and intermediate 
values of U, because the single particle band gap due to 
A is suppressed around edges, which enhances the inter- 
action effects. When U > U c « 4.9, the bulk antiferro- 
magnetism develops. The antiferromagnetic correlations 
along both the edge and central rows are enhanced by 
U. However, their difference is suppressed due to the 
disappearance of the helical edge states. 

The finite-size scaling of the xx-AFSF for the edge 
rows for different values of U are presented in Fig. [71 
Compared with the xx-AFSF calculated in the bulk 
(Fig. [2]), the edge antiferromagnetic correlations are 
much stronger than those of the bulk. Although the ex- 
trapolation to the infinite size in Fig. [71 implies a finite 
value of the Neel order of S x on the edge, we believe that 
it is an artifact due to the power-law scaling of the AF 
correlations. The ID nature of the edge states and the 
conservation of S z prohibits the true long range Neel or- 
dering of S Xj y but allows the quasi- long-range ordering, 
which is confirmed in the two-point spin correlations in 
Sect. HVCl 




• QMC data 




l=0A 11=1.0, ~ |r- 


T.a=1-05 


X=0A U=1.5, ~ |r- 


T,a=0.81 


a=0.1 U=2.0, ~ |r- 


T,a=0.64 



|r-r'| 



FIG. 8: (Color online) The two-point equal-time spin correla- 
tion functions along the zig-zag edge with A = 0.1 at values of 
U denoted in the insets. The sizes of the ribbon is 2 x 34 x 4. 
Because the zig-zag edge contains the sites of both A and B 
type, three different types of correlations are plotted in (a), 
(b), and (c), respectively. The Luttinger parameters are fitted 
from the correlation among A-sites on the tips as K w 0.8, 0.5 
and 0.4 for U = 1, 1.5 and 2, respectively. 



C. The stability of the helical edges 



According to the bosonization analysis in Ref. [l4j |, 
the scaling dimension of the 2k f Neel order of the xy- 
components is K, thus their equal- time correlations de- 
cays as l/\x — x'\ 2K . If the condition of the conservation 
of S z is released, a time-reversal invariant two-particle 
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correlated spin-flip backscattering term is allowed as 

Hbg,2 P ct = J dx ^j^dxip^ipLidxtpL^ + h.c. (17) 

At the particle-hole symmetric point of the KMH-model 
that we are simulating, the above term becomes the Umk- 
lapp term which conserves the lattice momentum. Such 
a term reduces the U(l) spin symmetry down to Z 2 . It 
has the scaling dimension 4if , and becomes relevant at 
K < K c = 1/2. In this case, it opens a gap by developing 
the long range 2k f magnetic ordering of S x or S y . Even 
for the cases that the two-particle spin-flip backscatter- 
ing are random disordered or at a single site, they still 
can destabilize the helical edge states at smaller values 
of the Luttinger parameter K [141 ]. 

According to the above analysis, the bulk paramag- 
netic regime at weak and intermediate coupling strengths 
should be divided into two regimes. At weak interactions, 
the helical edge states are stable against interaction ef- 
fects. The two-particle backscattering terms only have 
perturbative effects. On the other hand, at intermediate 
level of interaction strength, interaction effects are non- 
perturbative which breaks time-reversal symmetry along 
edges and thus destroys the helical edges. We emphasize 
that this destabilizing helical edges occurs when the bulk 
remains paramagnetic and time-reversal invariant. 

To numerically verify this picture, we present the cal- 
culation of the real space equal-time two-point correla- 
tions along the zig-zag edge in Fig. [8j Since each unit 
cell contains two non-equivalent sites, we denote the sites 
on the tips of the edge as A-sites and the other slightly 
inner sites as 5-sites. The correlation functions are de- 
fined as 

C AA (r,r') = (G\S^r)S^(f)\G), 
C BB (r,r') = (G\S*(r)S*(7)\G}, 

C AB (r,r') = 1 -{{G\S*{7)S*{7)\G) 

+ (G\S*(?)S£(7)\G)}, (18) 

where r and r* are along the zig-zag edge. The simulated 
results for A = 0.1 are plotted at different values of U 
in the bulk paramagnetic regime. The edge spin correla- 
tion exhibits the ferrimagnetic correlations among A and 
5-sites because the edge breaks the equivalence between 
A and 5-sites. The magnetic correlations are stronger 
among the outer A-sites, and are weaker among the in- 
ner B-sites. All of these correlations obey the power law 
and their decay exponents (a) are fitted. As further in- 
creasing U towards to the bulk antiferromagnetic regime, 
the difference between AA and BB correlations become 
weaker. 

Due to the domination of the magnetic correlation at 
A-sites, we use the decay exponents of Caa to fit the ef- 
fective Luttinger parameter K for the helical edge. The 
three plots in Fig. [8] (a) at U = 1, 1.5 and 2 gives rises 
to K = \a ~ 0.8,0.5, and 0.4, respectively. The case of 
U = 1 belongs to the topological band insulating phase in 




FIG. 9: (Color online) The definition of the positive direction 
for the NNN bonds in the honeycomb lattice, based on which 
the NNN current form factors Qq F and Qs F are denned in 
Eq. [H 



which interaction effects are perturbative. For the case 
of U = 2 at which the bulk remains non-magnetic, al- 
though the edge remains gapless, it is only because the 
conservation of S z which is not an essential symmetry 
of topological insulators. As long as the above Umklapp 
term Eq. [TT] is introduced, which unfortunately cannot 
be simulated by our QMC method, the gapless helical 
edge states are destabilized. We argue that the system 
enters a new phase with paramagnetic bulk but unstable 
edges. The transition point between these two paramag- 
netic phases at A = 0.1 lies at U ~ 1.5 with K « 0.5. 

We have calculated the edge spin correlations for other 
values of spin-orbit coupling and interaction parameters 
to map the boundary with K = 0.5 between two different 
bulk paramagnetic phases. The boundary is plotted in 
Fig. [3l As A decreases, the dispersion of the edge spectra 
becomes more flat, and interaction effects go stronger. As 
a result, the boundary shifts to lower values of U. In par- 
ticular at A = 0, the edge spectra become exactly flat, we 
expect edge ferromagnetism at infinitesimal U due to the 
density of state divergence. Thus the boundary should 
pass the origin. In particular, the edge ferromagnetism 
of graphene ribbon has been simulated in Ref. [38| . 



V. ABSENCE OF THE SPIN-ORBIT ORDER IN 
SPIN LIQUID PHASE AT A = 

Since Meng et al [36| claimed the existence of a spin- 
liquid phase for Hubbard model (X/t = 0) at 3.7 < U/t < 
4.3 (see Fig. [3j), it has attracted considerable interests 
and debates on the nature of this phase. One possibility 
of such a phase is that it could be a relative spin-orbit 
symmetry breaking phase with a non-trivial mean-field 
band structure [39[ . If it is the case, a finite X/t behaves 
like an external field to pin down the order parameter 
along the external spin-orbit configuration. Then the 
semi-metal and spin-liquid phase are indistinguishable at 
finite X/t. In this section, we will check the form factor 
of the such a spin-orbit order parameter between NNN 
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1/L 

FIG. 10: (Color online) The finite-size scaling of the form 
factors of the NNN currents in the charge sector Q^ F and 
Qs F in the spin sector as defined in Eq. [19] The periodical 
boundary condition is employed. The sample size is N = 
2 x L x L with L = 3, 6, 9 and 12. A is set in this calculation, 
and [7=3,4 and 5. 

sites at A = 0, and find negative results. 

Without loss of generality, we only consider the hori- 
zontal bonds. We define the positive directions for the 
NNN horizontal bonds as depicted in Fig. [9] Two dif- 
ferent NNN current orders are designed, including the 
charge flux order and the Kane-Mele type spin-orbit or- 
der, or, equivalently, the spin-current flux order. Their 
form factors are denoted as Q^ F and Qg F and are de- 
fined as 

Qc F = 7^|{B-l)^} 2 |G>, 

Qs F = ^(G\{J2(-^Ji,i + eS\G) (19) 

i 

where (— ) l takes the values of 1 or —1 for site i in 
the A and B sublattices, respectively; the charge cur- 

rent J ?,i+e x = J M^i!;t + J M+iL> and s P in current 
Ji,i + e x = ^™ ;t - J ™+; is the NNN vector along 
horizontal direction. Please note that the bond current 
operator here Jfi+g . a is different from that in Eq. [71 as 

J^l, <7 = i{c{ a c i+ ^-h.c.}, (20) 

where no summation over a is assumed. 

We have performed the simulation of the NNN charge 
and spin-current form factors defined in Eq. [19] for the 
Hubbard model at A = 0. The extrapolations of the 
form factors to the infinite lattice size are depicted in 
Fig. [10] The curves represents three typical Hubbard 
U values U = 3,4 and 5, which fall in semi-metal phase, 
spin-liquid phase and Mott insulating phase, respectively. 



For all the three parameters, both the charge and spin 
NNN current antiferromagnetic form factors vanish in 
thermodynamic limit, indicating the absence of the NNN 
charge and spin-current orders in all these three phases, 
especially the spin liquid phase. The nature of this spin- 
liquid phase, whether it is actually a subtly ordered phase 
or a genuinely exotic phase with non-trivial topological 
property, remains an unsolved question. 



VI. CONCLUSIONS 

We have studied the particle-hole symmetry in the 
KMH model, which results in the absence of the charge 
and spin currents and the absence of the quantum Monte- 
Carlo sign problem. The determinant QMC simulations 
have been performed for both the bulk and edge prop- 
erties. The bulk antiferromagnetic long range order ap- 
pears at large values of U. With the open boundary 
condition, the antiferromagnetic correlation is strongest 
along edges. 

We also studied the stability of helical edges in param- 
agnetic insulating phase when turn on infinitesimal two- 
particle backscattering term, which can be introduced 
by time-reversal invariant but S z not conserved interac- 
tion terms. The paramagnetic insulating phase in Fig. 
[3] can be classified into two regimes of weak and in- 
termediate interactions, respectively. In the weak in- 
teraction regime, the helical edge states remain gapless 
which is robust against the two-particle back-scattering; 
in the intermediate interaction regime, the edge states 
can spontaneously break time-reversal symmetry by de- 
veloping magnetic ordering along the edge by the two- 
particle backscattering term. Since this destabilizing he- 
lical edges occurs when the bulk remains time-reversal 
invariant, it is an interesting and open question whether 
the non-trivial bulk ^-topology is still maintained in this 
regime. 

We also checked that the spin-liquid phase in the Hub- 
bard model at A = in the honeycomb lattice is neither a 
spontaneously developed Haldane-type quantum anoma- 
lous Hall insulator, nor, the Kane-Mele type quantum 
spin Hall insulator. 
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Note Added During the preparation of this 
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manuscript, we learned the work on the QMC simula- 
tion on the same KMH model by Hohenadler et al [ioj . 



Appendix A: Absence of the sign problem for the 
zero temperature QMC simulations 

In this appendix, we prove the absence of the sign prob- 
lem of the KMH model at half-filling for zero tempera- 
ture determinant QMC, which is essentially an imaginary 
time projector algorithm. The explanation to the algo- 
rithm can be find in Ref. [37]. For readers' convenience, 
we also give a brief introduction below. 

The Hamiltonian composes of free and interaction 
parts 



H = H t + Hi. 

The free part reads 



(Al) 



(A2) 



where the kinetic energy matrix kernels and of 
the Kane-Mele model in Eq. IA2I are given in Section II. 
They satisfy the relation of Eq. [TTJ The interaction part 
is 



U 



where |^o) is the ground state; 9 is a projection param- 
eter large enough to ensure the trial wavefunction \iPt) 
is projected to the ground state |^o). The discretized 
HS transformation of the interaction term Eq. [A3] is per- 
formed in the density channel as the same as that in Eq. 
[9j The imaginary time propagator, i.e., the projection 
operator, is represented as 



6 -sh = ^2{u {l} (&,o)i[^ p (iy 



{i} 
l 



U {l} (e,0) = Y[ e" AT ^ °U K h°n e i V* i ¥X t ct t vtAi)cii 

1 



p=M 



(A5) 



where Ji, p (l) and rji iP (l) are the space-time discretized 
HS fields defined in Eq. [9] with I taking values of ±1, ±2; 

represents the summation over the spatial and tem- 
poral configurations of the HS field; /7{^(O,0) is the 
propagation operator for the HS configuration {/}. 

The trial wavefunction \i/jt) is required to be a Slater 
determinant, which we will specify later. The ground 
state |^o ) can be obtained from applying the imaginary 



(A3) time propagator e GH of Eq. IA5l on \i/jt) as 



The expectation value of a physical observable opera- 
tor O at zero temperature is defined as 

{6) = (MOM = (iP T \e- eH de-® H \ip T ) ^ (M) 



(V'olV'o) 



(^ T \e-^ H \4> T ) 



l^o> = E{ f/ m(©>0)n7 i , P (0e- l?? - ( ' ) ^ ¥ }lV'T), (A6) 
We further perform the calculation of Eq. IA4I as 



E w {(^t 1^(26,9) 6 UtieMMU^iiA^- 1 ^^} 
E w ^T|c/ ; (2e,o)|^ T )n,, p 7i )P (Oe- i " i - ( ' ) v^ 



J2P{i} (0) {l} , (A7) 

{»} 



r 



where (0){;} is the average value of O for the space-time 
HS configuration {1} defined as 

{0){l} - (^|c/ m ( 2 e,o)|^> ' (A8) 

and Pyy is the corresponding probability of the HS field 
configuration {/} as 

P{i} = \ (^1^(26,0)1^) n7i )P (0e- % ' p(0v ^- 

i,p 

(A9) 



Z is defined as Z = The summation over the 

HS configurations {/} can be done by using the Monte 
Carlo method. 



Next we prove the absence of the sign problem for the 
KMH model with purely imaginary NNN hoppings at 
half-filling in the zero temperature QMC method, i.e., 
the probability P{/} is positive-definite. We factorize 

the I^t) = ® IVv^)' wnere \^t^) i s a Slater- 

determinant state for spin-^ electrons with the particle 
number A/^, and similar convention applies for I^t^)- 
Then reads as 
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P{1} = f[ e -ArE i , 3 c l V<.c itei yA^E,*, P (0(cI t C T -|) |^t } 



p=2M 
1 



p=2M 



(A10) 



where the matrices if?- and K^- satisfy the relation of and 
Eq. [TTJ the HS fields Ji :P (l) are positive-definite. 

Let us perform a particle-hole transformation only to 
the spin-^ component 



c , 



^ = (-1)^, c 4 ^ = (-l)*c 4 , (All) 



then Slater-determinant state ) changes to another 
Slater-determinant state of holes with the hole number 
N-N^ denoted as ). We arrive at 



{1} 



p=2M 



x e i^ArU/2 V\ W,p (0(4t c «t " i ) 
X 



fc,iV-^| TT „-Ar Ei,j d li K ij d 3i 



n e- 

p=2M 



x e -t v /A7£V2E i *?i,p(0(diVu-i) ^Mr-JV^ 
x n^,p(0- (A12) 



det [Q t e M Q / ] , 



>1 "Mc 



n^^io) 



(A15) 



where M is an N x iV Hermitian matrix, or anti- 
Hermitian matrix. Based on these properties, we have 



{1} 



det 



x det 



(Qt) f ( II e- Kt e iVp(l) 

p=2M 
1 



Qt 



p=2M 



(A16) 



Now we add backfire explicit^ form of the Slater- where the matrix kernelg gatisfy Rt = {R iy and is 

a purely real diagonal matrix whose i-th diagonal element 



determinant states \ip^) ano ^ \*Pt' N as 



JV t N 



I ih,N-N l \ 



w) = n(E c M )J )io)=n(^ t ).i°)' 

■=i i=i j=i 3 

n (e 4^)10), 

= n (^).IO)^ (A13) 

3=1 3 

where |0) and |0)^, are the particle vacuum and hole vac- 
uum states, respectively; N is the number of of lattice 
sites; is a N x TV^-dimensional rectangular matrix, 
and is a N x (TV — N^)- dimensional matrix; and eft 
are vector notations for c\ and d| with i = 1 to iV. 

The Slater-determinant wavefunction has nice proper- 
ties as 



iV 



reads 



[*5>(0]< 



Ar ^- »?t,p(0- 



(A17) 



If we set the trial wavefunction to satisfy = = 
N/2 and = (Q T )*, then we have 



1 

Z 



det 



(Q T ) f ( n e- Kt e vp{i) ] ^ 

p=2M 



(A18) 



JJ^Q^-IO) = n^ eM< 3 t ]il°)' ( A14 ) at half-filling. 



thus the probability distribution is positive-definite 
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